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ABSTRACT 

We describe the iniplementation and testing of a smoothed particle hydrodynamics 
(SPH) code that solves the equations of radiation hydrodynamics in the flux-limited 
diffusion (FLD) approximation. The SPH equations of radiation hydrodynamics for an 
explicit integration scheme are derived and tested. We also discuss the implementation 
of an implicit numerical scheme for solving the radiation equations that allows the sys- 
tem to be evolved using timesteps much longer than the explicit radiation timestep. 
The code is tested on a variety of one-dimensional radiation hydrodynamics prob- 
lems including radiation propagating in an optically thin medium, optically thick and 
thin shocks, subcritical and supercritical radiating shocks, and a radiation dominated 
shock. Many of the tests were also performed by /Turner & Stone ( 200l|) to test their 
implementation of a FLD module for the ZEUS-2D code. The SPH code performs at 
least as well as the ZEUS-2D code in these tests. 

Key words: hydrodynamics - methods: numerical - radiative transfer. 



1 INTRODUCTION 

The smoothed particle hydrodynamics (SPH) technique was first introduced bv lLucvl ijlQTTT and lOingold fc Monaghanl ( 12Z3) 
as a Lagrangian method for solving astrophysical fluid dynamics problems. Over the subsequent years, it has been used to 
d ' investigate a wide range of astrono mical problems in cluding galaxy formation, star formation, stellar collisions, supernova 
explosions, and meteor impacts (see lMonagharJll992l for a review). It has also made its way into other fields of science and 
engineering, having been used to model tsunamis, volcanic eruptions, and ballistics. 

Many improvements and extensions to the method have been made over the years. Particle smoo thing lengths that 
varied in space as well as time were introduced iEvrardlllQSSl : iHernauist fc Katall989l : iBenz et al.lll99Cl) to allow the spa- 
tial resolution of the method to vary with density. This feature allows SPH to compete against the more modern adaptive 
mesh refinement (AMR) codes. Gravity, while naturally requiring a computational effo rt scaling as the square of the parti- 
cle n umber A'', was reduced to order A'' log A'' scaling using hierarchical tree structures IIHernauist fc Katdll989l : iBenz et alJ 
1Q20). Recently, a reformulation of the SPH equations has allowed bo th energy and entropy to be conserved simultaneously 

stl2n02l:lMonaghanl2nO!JI. Finally, many authors have suggested various reformu- 
UCingold fc Monaghanll983l : lEvrardll 988: Bal sara. Dll989l:lFlebbe et alll994l: 



ijinutsuka fc Imaedal20nit ISoringel fc H^rnq^ist|2(^(J2|;|M9n^ ghanl2nO!JI. Finally, many authors have suggested various reformu- 
lations of SPH's artificial viscosity iWoodll98ltlGingold fc Monaghanll983l : lEvrardll 988: Ba lsara. Dll989l:lFlebbe et alll994l: 
Watkins et a l. 1996VMor ris fc MonagharJl997^ or its elimination in favour of Godunov schemes Knutsukall994l : llnutsuka. S.- I.I 



19981: llnutsuka fc ImaedalboOlt ICha. S.- Hj|2002l:lcha fc Whitworthll200J) . 



However, with a few exceptions, SPH is primarily used to model pure hydrodynamics. Incorporating magnetic flelds into 
the SPH formalism has proven difficult. Problems encountered include numerical tensional instabilities, non-conservation of 
mome ntum, and/or not enforcing the divergence of the mag netic field to be zero (see th e introduction of IPrice fc Monaghanl 
|2004a|). Recently, a promising new method was presented bv lPrice fc Monaghanl ll2004allbr which may at last overcome these 
problems. ^_^^_^ ^^^_^^^^_^ ^_^ 

Similarly, SPH has been used for radiation hydrodynamics iLucylll977l : lBrookshawlll985l . 119861) . but not without encoun- 
tering problems. Boundary conditions that allow appropriate heat losses can be difficult to implement. More importantly, 
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the radiation equations involve the second spatial derivative of the radiation energy density. In the standard SPH formalism, 
the spatial derivative of a quantity is calculated using a sum over neighbouring particles, introducing noise. Calculating a 
second spatial derivative from such a first derivative produces a result that is very sensitive to particle disorder. In addition, 
for many interesting astrophysical problems, the timestep required to evolve the radiation equations explicitly is many orders 
of magnitude shorter than the dynamical timescale of the problem. Another recent approach to solving radiation transport 
with SPH is to combine SPH with Monte-Carlo radiative transfer ijOxlev fc WoolfsorJl200.'l) . This has the advantage that 
the radiative transfer is not carried out using the SPH formalism; the radiative transfer is calculated between hydrodynam- 
ical timesteps and the resulting temperatures are applied to the SPH particles. Since the method involves following photon 
packets, it is simple to include radiating point sources (e.g. stars) and scattering, and it allows frequency- dependent radiative 
transfer. Its main disadvantages are that it assumes the matter and radiation temperatures are the same and it is exceedingly 
computationa lly expensive, especia l ly in the optically thick regime. 

Recentlv. lClearv fc MonaghanI il999r published a method for modelling heat conduction with SPH. The heat conduction 
equation is similar in form to the equation for radiative diffusion. Cleary & Monaghan's method avoids the sensitivity to 
particle disorder by reformulating the heat conduction equation to use a first-order rather than sec ond-order spatial derivative. 
This s olves one of the hinderances for implementing radiative transfer within the SPH formalism. Ijubelgas. Springel fc Dolad 
l)2004) recently used this method to study heat conduction in galaxy clusters. 

In this paper, we describe the implementation of an SPH code that solves the equations of radiation hydrodynamics in 
the flux-limited diffusion approximation. The flux-limited diffusion approximation is briefly reviewed in Section 2. Section 3 
describes the SPH implementation in detail. Section 4 presents the results of various test calculations using the code. Our 
conclusions and some directions for future development are summarised in Section 5. 



2 EQUATIONS OF RADIATION HYDRODYNAMICS 

In a frame comoving with a radiating fluid, assuming local thermodynamic equilibrium (LTE), the coupled equations of 
radiation hydrodynamics (RHD) to order unity in v/c are 

■5f + /'V.« = 0, (1) 

p— ( — j = -V • F - Vt;:P + i-KKppB - cn^pE , (3) 

p— [-] = -PV ■ V - i-KUppB + CK.EpE , (4) 

pD(F\^^p_XpPp (5) 

c^ Dt \p J c ^ ' 

JMihalas fc Mihalaslll984l:lTurner fc StonelboOlh . In these equations, D/Dt = d/dt + w • V is the convective derivative. The 
symbols p, e,v, andp, represent the material mass density, energy density, velocity, and scalar isotropic pressure, respectively. 
The total frequency-integrated radiation energy density, momentum density or flux, and pressure tensor are represented by 
E, F, and P, respectively. These are the zeroth, first and second angular moments of the radiation specific intensity which, 
in general, varies with spatial position, viewing direction, frequency, and time. ^^^^^_^^^^^ ^_^ 

A full description of the fiux-limited diffusion approximation to the above equations is given by iTurner fc Stond (120011) . 
Here we simply summarise the main points. The assumption of LTE allows the rate of emission of radiation from the matter 
in equations 13 and 0] to be written as the Planck function, B. Equations |21 to |3 have been integrated over frequency, leading 
to the flux mean total opacity Xfj ^-nd the Planck mean and energy mean absorption opacities, kp and ke. In this paper, the 
opacities are assumed to be independent of frequency so that kp — ke and the subscripts may be omitted. The total opacity, 
X, is the sum of components due to abs orption, k, and scatter ing, a. Note that the above opacities have dimensions of length 
squared over mass (i.e. cm'^/g) whereas iTurner fc Stond 11200 ih deflne their opacities to have dimensions of inverse length (i.e. 
their opacities are equal to ours multiplied by the density). 

The equations of RHD may be closed by an equation of state specifying the gas pressure, the addition of consitutive 
relations for the Planck function and opacities, and an assumption about the relationship between the angular moments of 
the radiation fleld. In this paper, we use an ideal equation of state for the gas pressure p = (7 — l)up, where u — e/p is the 
speciflc energy of the gas. Thus, the temperature of the gas is Tg = (7 — l)pu/Rg^ = u/c^ where p is the dimensionless mean 
particle mass, Rg is the gas constant, and Cv is the specific heat capacity of the gas. The Planck function B — {a-B/'K)T^, 
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where ctb is the Stefan-Bohzmann constant. The radiation energy density also has an associated temperature Tr from the 
equation E — iaBT^/c. 

When the radiation field is isotropic P — ^E. The Eddington approximation assumes this relation holds everywhere and 
implies that, in a steady state, equation |3 becomes 

F = -^\/E. (6) 

This expression gives the correct flux in optically thick regions, where XP is large. However, in optically thin regions where 
XP — > the flux tends to inflnity whereas in reality |i^| ^ cE. Flux -limited diffusion solves t h is pro blem by limiting the flux 
in optically thin environment to always obey the above inequality. iLevermore fc Pomranind l|198lh wrote the radiation flux 
in the form of Pick's law of diffusion as 

F = -DV£, (7) 

with a diffusion constant given by 

D=!^. (8) 

XP 

The dimensionless function X{E) is called the flux limiter. The radiation pressure tensor may then be written in terms of the 
radiation energy density as 

P= fE, (9) 

where the components of the Eddington tensor, f, are given by 

f=i(l-/)I+i(3/-l)nn, (10) 

where n = Vi?/| Vi5| is the unit vector in the direction of the radiation energy density gradient and the dimensionless scalar 
function f{E) is called the Eddington factor. The flux limiter and the Eddington factor are related by 

f = X + X^R\ (11) 

where R is the dimensionless quantity R = \\7E\/{xpE). 

Equations 171 to ITTl close the equations of RHD, eliminating the need to solve equation |S] However, we must still choose 
an expression for the flux limiter, A. In this paper, we choose Levermore & Pomraning's flux limiter 

^(^) - 6tAtr^- ^''^ 

to allow direct comparison of our results with those of Turner & Stone (2001). In the optically thin limit 
hm X{R) = i (13) 

R — ^oo It 

to first order in R^^ so the magnitude of the flux approaches \F\ — c\\IE\/{xpR) = cE, as required. In the optically thick or 
diffusion limit 

lim X{R) = \ (14) 

to first order in R, so the flux takes the value given by equation |S| Many other forms of flux limiter have been developed to 
give more realistic performances on different problems (see Turner & Stone 2001 and references therein). 

3 NUMERICAL METHOD 

The radiation terms in equations|5]to|l]are added to a one-dimensional SPH code that is otherwise standard fe.g. lBenz et alJ 
ll99(l : lMonagharll992l ). The usual hydrodynamical SPH equations will not be rederived here. We use the standard cubic spline 
kernel, W , for the SPH summations over particles. The smoothing lengths of particles vary in time and space, subject to the 
constraint that the number of neighbours for each particle remains approximately constant at A'ncigh = 8. The smoothing 
length for particle i is given by 

hi = -MAX (l^i — a^i+sj 4- \xi — Xi+4,\, \xi — Xi-5\ + \xi — Xi-i\) , (15) 

where Xi is the position of particle i, and Xi+i and Xi+5 are the positions of the fourth and fifth closest particles in the postitive 
X direction, and Xi-4, and Xis are the positions of the fourth and fifth closest particles in the negative x direction. 

The SPH equations are integrated using one of two possible integrators: a predictor-corrector integrator, or a second- 
order Runge-Kutta-Fehlberg integrator. All particles are advanced with the same timestep. The integrators give essentially 
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identical results, so we only present the results of the test calculations usi ng one of the inteRrators, the second-order R unRe- 
Kutta-Fehlberg integrator. We use the standard form of artificial viscosity JGingold fc Monaghanlll983l : lMonaghanlll992l) with 
strength parameters a^ = 1 and /3v = 2 unless otherwise stated. The code does not include gravity or magnetic fields. 



3.1 Explicit SPH formulation of the radiation equations 

As SPH is a Lagrangian method, it is natural to use specific energies for the gas and radiation. The specific energy of the gas 
is w = e/p and we define the specific radiation energy to be ^ = E/p. Equation Q does not need to be solved directly since 
the density of each particle is calculated using the standard SPH summation over the particle and its neighbours. Equations 
121 to 12 can then be written 
Dv 



Dt P c 



Dt 
Du 



V ■ F Vv.V 



P 
pV • 



a 



P 



+ aCK, 



. (-) ) ' 

a Vcv/ 



(16) 
(17) 
(18) 



where a = 4(7b/c. The first terms on the right-hand sides of each of equations 1161 and 1181 are the hydrodynamic terms and are 
solved in the usual SPH manner (see below). 

The first term on the right-hand side of equation 1171 is the radiation flux term which, as discussed in Section |21 is given 
by 



VF 



P 



Kp 



P H \"-P / 

This has a similar form to the right-hand side of the equation for heat conduction 



dt p 



{kVT^ 



(19) 



(20) 



where k is the thermal conductivity. These equations both involve second-order spatial derivatives. As mentioned in the 
introduction, direct computation of a second-order derivative using standard SPH techn i ques makes the result extremely 
sensitive to particle disorder, and may result in unstable integration. IClearv fc MonaghanI (|l999|) overcome this difficulty by 
reformulating the second -order derivative as a first-order derivative using a Taylor series expansion. A full description of this 
reformulation is given bvl.Tu belgas et alJ l|2no4) . The heat conduction equation in SPH formalism then becomes 



'dT 



N 



4trC^ rC T 



pipj y rCi T" f^j 



vm, 



(21) 



where Wij = W{rij,h 

SPH formulation of equation 1191 is 



ij,;i,ij; with rij — Ti — Tj and hij = [hi -f hj)/2, and the subscripts denote particles i and j. Hence, the 



(£k 

\ Dt 



N 



j=i 



PiPj 



jPj 



jPj 



{piii - piij) 



VWr- 



(22) 



Because energy is transferred between pairs of particles, this formulation ensures that energy is conserved. 

The last terms in equations 1 1 71 and 1181 control the interaction between the radiation and the gas. Note that these terms 
are identical exce pt in sign, indicating that eri ergy is conserved. They are of the form T^ — Tg , as ^ = aT^ / p and u = CvTg, 
similar to that of iBlack fc Bodenheimeii l|197!Th except for our use of specific energies rather than energy densities. 

In one dimension, the radiation pressure term in equation 1181 becomes 

( — ) 

V -L^t- / radprcs 

where the Eddington factor is 

,2 / iv(M«)r ' 



--Vw 
P 



(V ■ v)^ M, 



(23) 



Ji — Ai + Xi 



i^ipiii 



The final SPH formulations of equations 1171 and 1181 are 

4- 






iPi K-jPi 



(^ 

V"'"' 



+ 



^0 



{Piii - Pi^i) 



VW^■ 




(24) 



(25) 



^4|:(|+|+''«)'"'«-'-»-»(^-(|-)') 
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(26) 



where Vij = Vi — Vj, and the first term in equation 1261 is the standard SPH expression for the hydrodynamical energy term 
p(V ■ v)/p when the thermodynamic variable of integration is energy. We note that the radiation terms can easily be expressed 
in terms of entropy if one wishes to add them to the SPH formalism of ISpringel fc Hernauisd 112002^ . We use the standard 
SPH artificial viscosity 



n„ = 



(— QvCs fiij + I3v fiij^) /pij if Vij- rij < , and 



if Vij ■ rij > 



where fiij — h [vi — Vj) ■ (r^ — rj) / (\ri — r_,| + 77^), with 77^ = O.Olft^ to prevent numerical divergences if particles get too 
close together. 

The full SPH expression of the momentum eauation ll6l is 



N 



j=i \ '■ 3 / j^i 



-^ = - ^ m J ^ + ^ + H,, VW{r,„ ft,,) - -1 ^ mj^.VWin,, h,,) , (27) 



Dt ^ ' \pj ' P^ ' 'V '■"'"' P, 

where the first term is the standard hydrodynamical term and the second term is due to radiation pressure. 



3.1.1 Explicit timestep criteria 

Explicit integration requires that the timesteps obey the Courant conditions for both the hydrodynamic and radiation pro- 
cesses. The usual hydrodynamical SPH timestep criteria are 

diCourant.i = \ \ 7 \ ^T , (28) 

Cs + fti |V • v\- + 1.2 (qvCs + P^hi |V • v\.) 
and 

diforce.z = Ca/ A-, (29) 

where we use a Courant number of (" = 0.3 and a^ is the acceleration of particle i. With radiation hydrodynamics, we modify 
equation 1281 so that the sound speed Cs used for computing the timestep includes the contribution from radiation pressure, 
and is given by Cs — [max (7, |j .&■"! 2 ^ where Ptot is the sum of both the gas and the magnitude of the radiation pressure, 
except where the sound speed is used as part of the viscosity, where it remains based on just the gas pressure. The radiation 
terms in the energy and momentum equations require additional constraints. The fiux diffusion timestep is 

di.ad,, = C%^- (30) 

C\i 

The timestep associated with the radiation pressure term in equation 1251 is simply 

diradprcs.i = C r „ , /v-7 ITT' ' ' 

fipie,i (V ■ 11), 

The timestep for the effect of radiation pressure on the momentum equation is taken care of already by the standard SPH 
force timestep criterion. The two timesteps for the interaction between radiation and matter are 

(32) 




(33) 



In cases where particles start off in thermal equilibrium (Tg — Tr), we found it was sometimes necessary to limit the 
timesteps further to ensure that too large a timestep was not used initially. In these cases, we also limited the timestep using 

, , , 4^ 

diint^.i 



CCi/ («CKi^), diint4,i = Cii/ ( acm ( ^ j I , dtintu.i = C^i/ (ficKi^) and dfmtu.i = C^i/ ( acm ( ^ j 
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3.2 Implicit integration 



Because the speed of light is typically much larger than the hydrodynamical velocity, the timestep obtained from the Courant 
condition for these radiative processes is much smaller than that for the pure hydrodynamics. It is therefore desirable to 
inte grate the radiativ e transfer implicitly so that a much larger timestep can be taken. 

iMonagharJ l|l997fl described an implicit integration scheme for handling the interaction between dust and gas particles in 
an SPH code. This is an iterative method that involves sweeping over all pairs of neighbouring particles a number of times. In 
Monaghan's case, the drag term between gas and dust particles was integrated implicitly. Here, the radiation and gas energy 
equati ons are integrated implicitly. Monaghan tried both backward Euler and two-step implicit integration methods lMonaghanI 
il997l section 4). We tried both of these methods and the standard trapezoidal method. We found the backward Euler and 
trapezoidal methods gave similar results, but that the two-step method did not converge in the limit of a large number o f 
sweeps. Since the trapezoidal method is the most accurate A-stable second-order linear multistep method llDahlauisdll963r . 
we chose to use this scheme. 

For the integration between a time t — n and t = n + 1 this scheme states that for a variable A{t), 



A. 



n-t-l 



,„ dt fdA? dA^+^ 



(34) 



where A" and A"^^ indicate the value of A at times t = n and t — n + 1, respectively. 

W e als o tried an iterative sche me based on the alternating-direction implicit (ADI) scheme used bv lBlack fc Bodenheimeil 
il975h and lTurner fc Stong i200J) . where progressively larger pseudo-timesteps are used between sweeps. However, testing 
showed that it was slower than the implicit method used below. 



3.2.1 Implicit flux diffusion 

Implicit integration of the radiation flux term can be performed by considering the interaction between two SPH particles i 
and j. According to the trapezoidal scheme 



^n + l _ f n . 1 dt 



— -b [piii - pjij + Piii 
piPj ^ 






(35) 



where b = 



for brevity, and N is the total number of sweeps. Similarly for particle j 



fn + l _ t'l I __ 






PiPj 



PiKi 



VW^„ 



It is possible to solve this set of equations for the quantity pif. 



n+l 



■ PjC^^. The solution is 



piii 



/-ri + l / ^Ti 

Pj?, = [Pi^i 



2 iv " nj \ Pi ^ pj ^ 



(36) 



(37) 



This is then substituted into equations 1351 and 1361 to provide the new values of the specific radiation energies of particles i 
and j. Performing a sweep involves calculating this interaction between all particle pairs. We tried both updating the energies 
after every interaction, and saving the changes from each interaction and updating the energies with them after all pairs had 
been considered. There was no significant difference between the two methods; we used the latter method since it minimizes 
truncation error. 



3.2.2 Implicit pressure and viscous force 



The pV ■ V and viscous contribution to the matter energy equation should also be done implicitly. Using the same technique 
as above, the interaction of particle j with particle i is 



n + l 



= Ui + 



Idt 

2iV 



1 f{i~l)K , (7-1)^" 



Pi 



Pi 



+ n„ ) mju.j VW,j + - ^ ^-^ + ^ -^ — + n„ ) m■jV,J\/W^J 

2 \ pi p-j 



,(38) 



and vice versa for the effect of i on j. Solving for u" /pi -\- it" / pj, we obtain 



,n+l 



+ ■ 



n + l 



Pj 



— + —] + 

Pi Pj 



(39) 



1^ ™v7w.. , (7-1) ^< , -n^ „,,.. , (7-1) /-: 



dt 
Npi 

dt 
Noting that VijVWij = VjiVWji, and that Hij = Ylji, this can be solved by defining four quantities 
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n + l 



U- 



2 4 \^ Pi Pj y 4: \ Pi Pj 

2 ^ \Pj pi J 4 \ Pj " 



'^r , di „ ^_ „„, , dt(7-l) ^< , u] 



Pi i/vpi N Api \ pi pj ' 



Pi ^Npj N 4pj \ Pi Pj 



c-i= — —^^ — 1njV^j\7W^j, 



dt (7 - 1) 
N Api 



dt (7 - 1) „,^^ 

C4 = — niiVrjVWij, 

N Apj 
and the solution for u^^^ / pi + u^^^ / pj is simply 

«, ^ /Pi + ttj /Pj = . (40) 

1 - C3 - C4 

This quantity is then substituted into equation 1381 and its counterpart for j to provide new values of Ui and Uj. Note that in 
fact the viscosity has a dependence on u, as the SPH viscosity involves the sound speed c^ ex y/u. However, we neglect this 
weak dependence during a sweep in order to simplify the mathematics, and simply update the sound speed in advance of the 
next sweep. We have not noticed any significant effect of this small approximation. 

3.2.3 Implicit radiation pressure 

The radiation pressure term in eauation l25l merelv involves the divergence of the velocity, and can easily be implicitly integrated 
thus: 

cr+' = er - i^ ( V • «), (cr+' + &) /., (41) 

whose solution is 

^' "^n(i + ^f(v-«),/0| ^^ 

3.2.^ Implicit radiation-matter interaction terms 

The implicit solution of the interaction terms between ^ and u for a given particle i involves the solution of a quartic 
(fourth-order polynomial) equation, which is non-trivial. The equations for ^ and u are 






/-n + l /-n 1 QC / Pil^i ^ri ^i / n\4 , Pi^i j-n-\-\ ■ -"i / n-|-i\^ i /^o\ 

C. =^» "2iv"'^i^^" ~^^"'^ +~^^' -I^("« ) ) (43) 



and 

n + l _ n I 1 dt ( pjKj ^„ Ki , „>4 Pi>^^n + 1 J<d ( n + l\4 



„ ± Ul, / pil-vi „ rvi „ 4 Pii^i ^n + 1 l^i ( n + l\*l /..n 



Here the term to solve for is 

7i+l\4 



^ = PiCi ^ ("» j ■ (45) 

Defining the quantities 



Cl = piii - 2jj'°-'^P^ [ ~^^^ ~ — (^i) 



Cl 



Hi a 



Cva'^^ 



C3 = Ui + -^jj^ac I — Ci - 7:^ (^^i ) 



i Ql / pitli ^n ^i / n\4 
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1 di Hii 

1 df Hi 

2 7!""^' 



the equation to be solved becomes 



(46) 



X^+'-^X-^+'4x'+('-^ + ^ + ^]x+^^i^i^=0, (47) 

which can be solved using any method for solving quartic equations (see Appendix^ for the analytical method). While it is 
possible to use a numeric solver, e.g. a Newton-Raphson method, we have found that the analytic solution is to be preferred 
since tests show it is much quicker. One minor difficulty when calculating the analytic solution is that while two of the roots 
of a quartic equation (usually) contain an imaginary component and can readily be discarded, the code must be able to 
distinguish between the two real roots. We choose the correct solution as the one which gives a new gas temperature that lies 
between the old values of Tg and Tr. In the case of four real roots, the same test can be applied. Four complex roots usually 
indicates that the timestep is too large. 

If the real root chosen is denoted by Xi, the quantity ^X\ should be substituted back into equations 1431 and 1441 as 
the term ( £i£i^^+i Ji- (u"+^ j J to provide the new radiation and gas energies. We have found that solving for X\ as 

piC"^^ — ■;f--x- y'^^) rather than I ^pii^^^ J^ (u"^^) I can be more accurate if the numbers involved in the calculation 

are closer to the order of unity, ensuring that the range does not fall outside the Fortran double precision variable bounds. 
For different values of opacity it may even be desirable to solve for ^X or even (^) X depending on the problem. The code 
can be made adaptive in this respect by checking that the co-efficient of the first-order term in equation 1471 is smaller than 
the fourth root of the largest allowed floating point number, and adapting the quantity to solve for accordingly. 

The analytical solution of this quartic equation, whilst being much faster than a numerical solver, can fail when the 
number of sweeps is too small. In this case, the code detects a failure to reach a physical temperature and returns to the 
calling subroutine requesting a greater number of sweeps. This is a fairly common occurance when the code is just starting, 
as the initial number of sweeps is an estimate input by the user. 



3.2.5 The Adaptive Sweeping Scheme 

The way the sweeping is performed is important, and can greatly affect the accuracy of the implicit method. Initially the 
quantities ^ and u are copied into arrays upon which the sweeping operations are to be performed. Then sweeping commences 
in the following order: First the flux-limiter A^ is calculated for each particle, along with the Eddington factor fi. The next 
step is to calculate the flux diffusion term for all i — j particle pairs, update ^, then calculate the radiation pressure term 
and update ^ again. Independent of the radiation energy calculations, the sound speeds Cs (used in the SPH viscosity) are 
calculated, followed by the gas pressure/viscosity term and u is updated. The sweep ends with the interaction term between 
u and ^, after which the next sweep commences with the calculation of the flux-limiter using the new values of ^ and u. Once 
the specified number of sweeps are complete the temporary arrays are copied back into their normal counterparts. 

Since it is not known a priori how many sweeps are required to reach a certain accuracy, we use an Adaptive Sweeping 
Scheme (ASS) to determine the necessary number of sweeps. The code begins by trying an initial number of sweeps 2™ (where 
m is an integer, m > 0). The implicit subroutine is called twice, once for 2™"^ sweeps and once for 2™ sweeps. If the fractional 
error in both ^ and u for all particles between the two sweeps is less than a specified tolerance the values of £, and u from the 
2™ sweep calculation are used, and m is decreased by unity for the next timestep to allow the code to adapt if fewer sweeps 
are required. We found that a tolerance of 10~ gave a good trade off between accuracy and speed and this value is used for 
all the test calculations presented in this paper. If the error is above the tolerance, m is incremented by unity, and the method 
repeats until either an acceptable solution is found, or too many sweeps are attempted. The trapezoidal scheme is such that 
doubling the number of sweeps typically halves the errors. The ASS can be optimized to take this into account, e.g. if the 
maximum error is one quarter of the acceptable tolerance the number of sweeps can be reduced by a factor of four for the 
next timestep. If at any time the implicit integration would take ^ or w negative, the implicit routine immeditately returns to 
ASS with a message to double the number of sweeps, similar to what occurs when the quartic solver fails. 

It is also desirable for two-step integrators like the Runge-Kutta-Fehlberg or predictor-corrector integrators that require 
the calculation of quantities at both the half and the full timesteps to use two independent numbers of sweeps. Typically, the 
full timestep call of the implicit method requires twice as many sweeps. 
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Figure 1. The evolution of the gas energy density e as it equiUbriates with a radiation energy density E = 10 erg cni~ . In the upper 
case, e = 10^" erg cm~^ initially, while in the lower case e = 10^ erg cm""^. The solid line is the analytic solution, the crosses are the 
results of the SPH code using implicit timesteps of the greater of lO"^"* s or five percent of the elapsed time, and the squares with a 
timestep of the larger of 10^^^ s or five percent of the elapsed time. The symbols are plotted every ten timesteps. 

4 TEST CALCULATIONS 

4.1 Heating and cooling terms 

As in lTurner fc Stong ||2001j) we tested the interaction between the radiation and the gas to check that the temperatures of 
the gas and the radiation equalise when Tg 7^ Tr initially. A gas was set up so that there was no initial velocity, with a density 



10 ^ g cm ^, opacity k 



0.4 cm g , and 7 



and E 



10^^ ergs cm ^. 



Two tests were carried out, one where the 

10^ ergs cm~'^. 



gas heated until it reached the radiation temperature, and one where it cooled. The first test had e = up 

and the second e = up — 10^" ergs cm""^. The boundaries of the calculation used refiective ghost particles, and the position 

of the boundary was fixed. 

This problem can be approximated by the differential equation 



At 



pc^ 



(48) 



in the case where the energy in the radiation is much greater than that in the gas. In figure Q the solid line is this analytic 
solution, plotted both for the cases where Tg increases and decreases. The crosses are the results of the SPH code using an 
implicit timestep that is set to the larger of 10"^** s or five percent of the time elapsed. The squares are similar, but with 
a timestep being the greater of 10^^^ s or five percent of the time elapsed. As can be seen, the match between the analytic 
solution and the solutions given by the SPH code is excellent. 



4.2 Propagating radiation in optically thin media 

In the standard diffusion approximation, radiation can, in optically thin regions, approach infinite velocity. This is unphysical, 
and the flux limiter has been introduced to limit the diffusion of radiation to the speed of light. To examine how well our code 
limits the speed of the radiation, a one centimetre long one-dimensional box is filled with 100 equally spaced SPH particles, 
with E = 10~^ erg cm~^ (^ = 0.4 erg g~^), p = 0.025 g cm~^ and k = 0.4 cm^ g~^. Initially, the radiation and gas are in 
thermal equilibrium. 

At the start of the simulation, the radiation energy density for the leftmost ten particles was changed to -E = 0.1 erg 
cm^'^(5 = 4 erg g~^) and the resulting radiation front was allowed to propagate across the region. The ghost particles were 
reflective except in speciflc radiation energy ^, which was flxed equal to ^ = 4 erg g~^ at the left hand boundary and S, = 0.4 
erg g~^ at the right. The implicit code was used with various timesteps ranging from an explicit timestep to a single implicit 
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Figure 2. The propagation of a radiation pulse across a uniform medium. The time is t = 10~ s. The vertical dashed line shows the 
expected position of the pulse based on the speed of light. The results are essentially independent of the size of the implicit timestep used. 
Results are given for implicit timesteps equal to (solid line), ten times (dotted line), and one hundred times (dashed line) the explicit 
timestep. The dot-dashed line gives the result using a single implicit step of 10~^^ s. The results for the explicit code are identical to 
that of the implicit code using the explicit timestep, and hence are not plotted. The initial conditions are also shown (solid line). 



step lasting 10~^^ s. The results are shown in figure|5| The result for the explicit code lies on top of the result for the implicit 
code using an explicit timestep. 

As can be seen, the radiation pulse propagates at the correct speed, even using a single implic it timestep that i s mor e 
than 10^ times the explicit timestep. The front is smoothe d out in a manner simil ar to the results of iTurner fc Stond 11200 If) : 
both methods are quite diffusive in this situation. We note lTurner fc Stond (120011) began with an initial discontinuity in the 
radiation energy density of 21 orders of magnitude. We found that as the discontinuity increases in magnitude, the radiation 
pulse propagates more slowly with SPH due smoothing of the gradient of E and the consequent error in the value of the 
flux-limiter. 



4.3 Optically-thick (adiabatic) and optically-thin (isothermal) shocks 

A shock-tube test was set up to investigate the way the code simulated optically-thin and optically-thick regimes and the 
transition between them. In the limit of high optical depth, the gas cannot cool because the radiation is trapped within the 
gas; thus the shock is adiabatic. An optically-thin shock, on the other hand, is able to efficiently radiate away the thermal 
energy and thus behaves as an isothermal shock. In these shock tests, the gas and radiation are highly coupled and, thus, 
their temperatures are equal. 

A domain 2 x lO'^^ cm long extending from x = — 1 x 10^^ to 1 x 10^^ cm was set up, with an initial density of p = 10"^" 
g cm""^, and the temperatures of the gas and radiation were initially 1500 K. One hundred particles were equally spaced in 
the domain, with those with negative x having a velocity equal to the adiabatic (7 = 5/3) sound speed vo = Cs = 3.2 x 10^ cm 
s~^, and those with positive x travelling at the same speed in the opposite direction. The two flows impact at the origin, and 
a shock forms. Opacities of «; = 40, 0.4, 4.0 x 10~^ and 4.0 x 10~^ cm^ g~^ were used to follow the transition from adiabatic to 
isothermal behaviour. Ghost particles were placed outside the boundaries and maintain the initial quantities of their respective 
real particles. The boundaries moved inwards with the same velocities a s the initial velocities of th e two streams. 

The adiabatic and isothermal limits can be solved analytically (e.g. lZel'dovich fc Raizerll2002f) . The shock speed is given 
by 



D = 



(7off - 3) + ^(7cff + ly^vl + 167e. 



(49) 
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Figure 3. A set of shocks with differing opacity at time t = 1.0 X lO' s. Density is on the left, and gas temperature on the right. The 
crosses are the SPH results; the solid line gives the analytic solution for an adiabatic shock, and the dashed line for an isothermal shock. 
The opacities are (top) k = 40, (upper middle) k = 0.4, (lower middle) ft = 4.0 X 10""^ and (bottom) k = 4.0 X 10~^ cm^ g"'^- As the 
opacity is decreased, the shocks transition from adiabatic to isothermal behaviour. 



where 7cff = 5/3 for the adiabatic case, and 7eff = 1 for the isothermal case. The ratio of the final to the initial density is 
given by 



(50) 



(51) 



po D 

and, for the adiabatic shock, the ratio of the final to the initial temperature is given by 

]j_ _ P]_ , ""pPpo 
To po Po 

These analytic solutions are shown by the solid and dashed lines in figure |3 In the figure, the opacity decreases from top to 
bottom showing the transition from optically-thick (adiabatic) to optically-thin (isothermal) behaviour. The extremes are in 
good agreement with their respective analytic adiabatic and isothermal solutions. Note that the spike in thermal energy near 
the origin and the corresponding reduction in density for the optically-thick case (due to 'wall-heating') is softened by the 
radiation transport that occurs in the intermediate opacity calculation with k = 0.4. 

Where possible, comparison shows that the results from a fully explicit calculation are virtually identical to those obtained 
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Figure 4. The sub-critical siioclc, witii piston velocity 6 x 10^ cm s~^ and 100 particles. Tlie large panels show the results using 
the explicit code. The top panels show radiation (solid line) and gas (dotted line) temperatures. The bottom left panel shows the 
normalised flux and bottom right panel the Eddington factor. The long-dashed lines give the analytic solutions for the gas temperature 
and normalised flux. An Eddington factor of 1/3 is also indicated for reference (short-dashed line, lower right panel). The subpanels plot 
the logarithm of the difference between the results using the explicit code and the implicit code (see the main text). The implicit code was 
run with timesteps of 1 (solid line), 10 (short-dashed line), and 100 (dotted line) times the explicit timestep and with a hydrodynamical 
timestep (long dashes). 

implicitly. However, as the opacity is increased, the implicit code offers an enormous improvement in speed because the explicit 
timestep is set by the interaction term between the radiation and gas (equations 1321 and I33II which is inversely proportional 
to the opacity. With an opacity oi k = 4.0 x 10~^, the codes are similar in speed. However, the implicit code is ~ 100 times 
faster than the explicit when k = 4.0 x 10"'^, and ~ 10* times faster when k = 0.4. For k = 40, we could not run an explicit 
calculation, but the implicit code is likely to be ~ 10® times faster. 



4.4 Sub- and super-critical shocks 

A supercritical shock occurs when the photons generated by a shock have sufficient energy to preheat the material upstream. 
The characteristic temperature profile of a supercritical shock is where the temperature on either side of the shock is sim- 



SPH with flux-limited diffusion 13 



^ 4000 




^ 4000 



M -4 

o 2x10'° 4x1010 6x1010 



Position, X / cm 




o 



10 -5 

Optical depth, t 




10 -5 

Optical depth, t 





1 


1 1 1 1 1 1 


1 1 1 1 1 1 


«4-. 








C 


0.8 


- 


— 


o 








-tJ 






~ 


o 




- 


- 


Id 


6 


— 


— 


u. 






- 


C 


7^ A 


o 
-1-1 


4 


- \ 





CU) 






"^ — 


C 


- 


- 


Tl 




- 


- 


T3 


OiH 


— 


— 


Ixi 




- 


- 




n 


1 1 1 1 1 1 


1 1 1 1 1 1 




10 -5 

Optical depth, t 



Figure 5. The super-critical shock, with piston velocity 1.6 X 10® cm s ^ and 100 particles. This shock is strong enough for radiation 
from the shock to preheat the gas upstream. See figure|4]for more details. 



ilar, rather than the dow nstream temperature being much higher than the upstream, as occurs in a subcritical shock (see 
IZel'dovich fc Raizeill2002l . for more details). ^^^^^^^^^_^^^^^^^^^ ^_^ ^^^^^_^^^^^ ^_^ 

The initial conditions of this problem are those of ISincell. Gehmevr fc Mihalaa 1119991) and lXurner fc Stonj 11200 If) . A gas 
with opacity k = 0.4 cm^ g~^, uniform density p = 7.78 x 10"^" g cm~^, mean molecular weight fi = 0.5 and 7 = | is 
set up with ^ and u in equilibrium, with a temperature gradient of T = 10 + (75a;/7 x 10^") K. Initially the particles are 
equally spaced between x = and a; = 7 x 10^" cm for the supercritical shock, and between a; = and x = 3.5 x 10^" cm 
for the subcritical shock. At time i = a piston starts to move into the fluid from the left-hand boundary (simulated by 
moving the location of the b oundary). For the su bcritical shock the piston velocity is tip = 6 km s~^, and for the supercritical 
shock Vp = 16km s~i, as per lSincell et al.l (|l999|). The ghost particles are reflective in the frame of reference of the boundary. 
Artificial viscosity parameters 0^ — 2 and /3v = 4 were used to smooth out oscillations. 

The results of calculations using 100 particles are shown in figure |1] for the sub-critical shock and in figure |3 for the 
super-critical shock. The top left panel for each shows the temperature of the radiation field (solid line) and the gas (dotted 
line) against position, and the top right shows the same quantities against optical depth r, with t — set at the shock front 
(measured from the density distribution). The botto m left panel shows n orm alised fiux, and the bottom right the value of the 
Eddington factor. The analytic solutions discussed bv lSincell et alJ 111999^ and lZel'dovich fc Raizerl (120021) for the temperatures 
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Figure 6. The sub-critical shock with differing resolutions. For 50 particles, the dot-long-dashed lines give the radiation temperature, 
the flux, and the Eddington factor while the short-dash-long-dashed lines give the gas temperature. For 100 particles, the solid lines 
give the radiation temperature, the flux, and the Eddington factor while the dotted lines give the gas temperature. For 250 particles, 
the dot-short-dashed lines give the radiation temperature, the flux, and the Eddington factor while the short-dashed lines give the gas 
temperature. The analytic solution is given by the long-dashed lines. 



and fluxes of the shocks are shown with long-dashed hnes. Figures 0] and |S] are plotted using the explicit code. In subpanels 
beneath the main panels, we compare the results from the implicit code with the explicit results. Calculations were performed 
using the implicit code with timesteps of 1, 10, and 100 times the explicit timestep and with a timestep limited only by the 
hydrodynamical timestep criteria (equations 1281 and 1291 . In the subpanels, we plot the differences of the implicit results with 
respect to the explicit results. We divide the difference between the implicit and explicit values by the explicit value to obtain 
a fractional error and take the logarithm of the absolute value of this fraction. Thus, a difference of —2 on the subpanels 
corresponds to an error of 1 percent with respect to the explicit result. 

Figures|H]E|and|5]show the effect on the shocks of changing the resolutions. For the sub-critical shock, the results converge 
towards the analytic solutions for both the gas temperature and the flux as the particle number is increased. Note that our 
highest resolution case has the same number of particles per unit length as lTurner fc Stond 11200 If) have grid cells. Increasing 
the resolution of the super-critical shock has less of an effect on the results, although the flux in the vicinity of the shock 
improves, and the spike in gas temperature at the location of the shock is better resolved (c.f. figure|HJ. 

Both the explicit and implicit codes model the sub-crit i cal an d super-critical shocks well. With similar reso lution, our 
results are in good agreement with those of iTurner fc Stond (|2001j), although we note that lTurner fc Stond (|2001j) only used 
timesteps as large as ten radiation diffusion times. However, we found that in these tests, our implicit code offers little 
advantage over the explicit code. With the super-critical shock and an implicit timestep 10 times longer than the explicit 
timestep, the speed is similar, but using a hydrodynamic timestep the code is roughly 10 times slower. For the sub-critical 
shock, the situation is worse. Running with an implicit timestep 100 times longer than the explicit timestep takes ~ 60 times 
longer and implicit calculation using a hydrodynamical timestep takes « 400 times longer than an explicit calculation. The 
decrease in speed is due to the large number of sweeps required by the implicit integration scheme in order to achieve the 
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Figure 7. The super-critical slioclc witli diflfering resolutions. For 100 particles, the dot-long-dashed lines give the radiation temperature, 
the flux, and the Eddington factor while the short-dash-long-dashed lines give the gas temperature. For 200 particles, the solid lines 
give the radiation temperature, the flux, and the Eddington factor while the dotted lines give the gas temperature. For 500 particles, 
the dot-short-dashed lines give the radiation temperature, the flux, and the Eddington factor while the short-dashed lines give the gas 
temperature. The analytic solution is given by the long-dashed lines. 
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Figure 8. A zoom-in of the spike for the supercritical shock. We plot the results with 100 (solid line), 200 (dotted line) and 500 (dashed 
line) particles. As expected, the spike is more obvious with higher resolution. 
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Figure 9. The radiation-dominated shock at time t = 5 X 10~^ s. We plot the velocity, radiation and gas energy densities and gas 
density versus position. The vertical dashed lines show the expected shock thickness. Gas flows into the shock from the left. The after 
effects of the transient that occurs at the start of the calculation can be seen at the far right of the plots. 



required accuracy. Increasing the tolerance above 10 "^ results in quicker calculations, roughly inversely proportional to the 
tolerance, however the solution is less accurate and both the temperatures of the gas and radiation are overestimated. 



4.5 Radiation-dominated shock 

In material of high optical depth the radiation generated in a shock cannot diffu se away at a high rate, and so the radiation 
becomes confined in a thin region adjacent to the shock. iTurner fc Stond (120011) performed a calculation that tests whether 
the shock thickness is what one would expect in these circumstances. An extremely high Mach number shock (Mach number 
of 658) is set up, with the gas on the left having an initial density of p = 0.01 g cm~"^ (set by using different mass particles), 
opacity k — 0.4 cm^ g~^, temperature T^ = Tg = 10* K, and speed 10^ cm s~i. The gas on the right has density p — 0.0685847 
g cm~^, opacity k — 0.4 cm^ g~^ , temperature T^ — T^ = 4.239 x 10^ K, and speed 1.458 x 10* cm s~^. The locations of 
the boundaries move with the same speed as their respective particles, and the properties of the ghost particles outside these 
boundaries are fixed at their initial values. We use the implicit code with a timestep 100 times larger than the exphcit timestep. 
1500 particles are equally spaced over a domain extending from x = —6 x 10^ cm to a; = 1.5 x 10^ cm initially, with the 
discontinuity aX x — 0.5 x 10 cm. The location of the shock should be fixed in this frame, although individual particles flow 
through the shock. 

After a period where a transient feature forms at the shock front and drifts downstream with the flow, a stable shock 
is established. Its thickness is expected to be roughly equal to the distance I = cX/npui, where Wi is the speed of material 
flowing into the shock front. 

Figure El shows the results at a time i = 5 x 10~ s, of the radiation energy density E, gas energy density e, velocity 
V, and density p versus position. The vertical dashed lines indicate the expected shock thickness and the SPH results are in 
good agreement. The after-effects of the transient moving downstream can be seen on the right of the plots of density and 
gas energy density. The explicit code takes roughly a factor of ten times longer than the implicit code for this test. 
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5 CONCLUSIONS 



In this paper, we have presented a method for including radiative transfer in the flux-limited diffusion approximation in a 
smoothed particle hydrodynamics (SPH) code. The energy equations may be integrated either explicitly or implicitly. By 
integrating the energy equations implicitly the short timesteps due to the rapid propagation of radiation can be avoided and 
timesteps limited only by hydrodynamical processes can be used. In problems where the temperatures of the gas and radiation 
are equal, the implicit code loses little in accuracy but may be many orders of magnitude faster than an explicit calculation. 
However, in cases where the two temperatures are not equal the implicit scheme may actually be slower due to the large 
number of sweeps required to obtain an accurate solution. 

Overall, the SPH code per forms as accurat e ly as the ZEUS code with radiative transfer in the flux-limited diffusion 
approximation as described by iTurner fc Stond 11200 ih . Although the code and test calculations presented here are one- 
dimensional, it is straightforward to implement the radiation hydrodynamical equations in a three-dimensional SPH code due 
to the fact that the SPH method is structured around interactions between pairs of particles. In going from one to two or 
three dimensions, the general method is the same. The main difference is that the radiation pressure term in the evolution 
equation for the specific radiation energy involves a more complicated tensor equation. Aside from this, the major limitation 
associated with increasing the number of dimensions is that more particles are required to obtain the same resolution, thus 
more computational effort is required. Three-dimensional calculations will be presented in a subsequent paper. We plan to use 
this method to investigate the effects of radiative transfer on star formation. Other applications include modelling phenomena 
such as accretion discs and radiation-dominated flows. 
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APPENDIX A: THE ANALYTIC SOLUTION OF A QUARTIC EQUATION 

To solve the quartic equation x'^ + a^x"^ + a-zx^ + aix + uq — Q, it is necessary to first find the real root j/i of the resolvent 
cubic 

y — a2y + (aiaa — 4ao)y + (4a2ao — fli — Ogao) = , (Al) 

either analytically or numerically. The algorithm for a real root (of which there must be at least one) can be obtained from a 
mathematical manipulation package and output as computer code to be copied into a program. The four roots of the quartic 
are then given by the four solutions of 



2 I 

Z + 



0.3 ^ f al 
yT^+yi-a2 



z + 



m 



y =F ( — I - ao 



: , (A2) 



JAbramowitz fc Stegunlll973) . If some of the coefficients are large and of roughly equal size it may be necessary to Taylor 
expand the relevant portion of the equation to maintain accuracy. 
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